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Abstract. 

Dark matter direct and indirect detection signals depend crucially on the dark 
matter distribution. While the formation of large scale structure is independent of the 
nature of the cold dark matter (CDM), the fate of inhomogeneities on sub-galactic 
scales, and hence the present day CDM distribution on these scales, depends on 
the micro-physics of the CDM particles. We study the density contrast of Weakly 
Interacting Massive Particles (WIMPs) on sub-galactic scales. We calculate the 
damping of the primordial power spectrum due to collisional damping and free- 
streaming of WIMPy CDM and show that free-streaming leads to a CDM power 
spectrum with a sharp cut-off at about 10“®Mq. We also calculate the transfer function 
for the growth of the inhomogeneities in the linear regime, taking into account the 
suppression in the growth of the CDM density contrast after matter-radiation equality 
due to baryons and show that our analytic results are in good agreement with numerical 
calculations. Combining the transfer function with the damping of the primordial 
fluctuations we produce a WMAP normalized primordial CDM power spectrum, which 
can serve as an input for high resolution CDM simulations. We find that the smallest 
inhomogeneities typically have co-moving radius of about 1 pc and enter the non-linear 
regime at a redshift of 60±20. We study the effect of scale dependence of the primordial 
power spectrum on these numbers and also use the spherical collapse model to make 
simple estimates of the properties of the first generation of WIMP halos to form. We 
And that the very first WIMPy halos may have a significant impact on indirect dark 
matter searches. 
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1. Introduction 

Analysis of the anisotropies in the cosmic microwave background (CMB) radiation P] 
hnds that the relative matter density = 0.29 ± 0.07 is signihcantly larger than 
the relative baryon density fib = 0.047 ± 0.006. This is consistent with the observed 
abundances of light elements and primordial nucleosynthesis (see e.g. Reference j2]) 
and the power spectrum found from galaxy red-shift surveys jS], and indicates that 
the Universe contains a signihcant amount of non-baryonic cold dark matter (CDM). 
The identihcation of the nature of the CDM particles is a major outstanding issue 
for cosmology (and would provide reassuring conhrmation of the CDM cosmological 
paradigm). 

There are various CDM candidates (for a recent review see jl]), the most well 
studied of which are Weakly Interacting Massive Particles (WIMPs) and axions. WIMPs 
are a particularly attractive CDM candidate, since a stable relic from the electroweak 
scale generically has an interesting present day density, flwimp ~ 0(1) 0 . There are a 
large number of ongoing experiments attempting to detect WIMPs directly in the lab jB] 
or indirectly via their annihilation products (gamma-rays, antiprotons and neutrinos) 
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The signals expected in dark matter detection experiments depend in most cases 
on the distribution of the dark matter. Direct detection experiments probe the dark 
matter distribution on sub-milli-pc scales IHli while indirect detection signals are 
strongest from the highest density regions of the Milky Way (see e.g. References nm) 
and the extra-galactic gamma-ray signal depends on the dumpiness of the WIMP 
distribution HH. Reliable predictions for the expected signals therefore require an 
understanding of the dumpiness of dark matter on small (sub-galactic) scales. The 
density perturbations on very small scales, and hence the properties of the hrst 
generation of structures to form, depend on the microphysics of the CDM and the 
present day density distribution may retain traces of these hrst structures. 

Two of the present authors showed that the microphysics of WIMPs (specihcally 
collisional damping due to interactions with the radiation component and free- 
streaming) lead to a fundamental small scale cut-off in the WIMP density perturbation 
power spectrum [la cni. Subsequently we presented the small scale WIMP density, 
velocity and potential perturbations and estimated the properties of the hrst generation 
of WIMP halos to form for the case of the WIMP being a bino M- Meanwhile 
Berenzinsky et al. ra studied the survival probability of these hrst halos analytically. 
More recently Diemand et al. jlH] have carried out high resolution simulations of the hrst 
WIMP halos to form and there has been (inconclusive) discussion about whether these 
halos will be disrupted by interactions with stars na. Shortly after the hrst version 
of the present work became available on arXiv.org, Loeb and Zaldarriaga published a 
numerical calculation of the cut-oh scale Their work is an improvement on the 

analytic treatment in the present paper, although the estimates for the cut-oh mass 
scales agree up to a factor of order unity (the diherence between their and our 
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1 O“®M 0 stems largely from different assumptions for the kinetic decoupling temperature 
- 10 MeV in versus 30 MeV for our benchmark model in ^3], which is a somewhat 
more realistic value for the lightest neutralino). 

In this paper we present the detailed calculations behind the results presented in 
our earlier letter m generalised to generic WIMPs, including a new more intuitive 
collisional damping calculation using the Navier-Stokes equation. We also compare our 
analytic transfer function with numerical calculations and study the effects of a scale 
dependent primordial power spectrum. In Section |21 we estimate the temperature scales 
related to chemical and kinetic decoupling for generic WIMPs. In Section El we dehne 
the fluid perturbation variables and present the equations necessary for our subsequent 
calculations and solve the evolution of the matter and radiation perturbations during 
the radiation dominated regime. In Sections 0] and El we calculate the damping of 
the WIMP density contrast due, respectively, to collisional damping (as a result of 
elastic interactions with other species) and free-streaming. We then (Sectional) calculate 
the transfer functions, which encode the evolution of the perturbations, by solving the 
evolution equations around matter-radiation equality and matching these solutions to 
the sub-horizon limits of the radiation domination solutions found in Section El and 
compare our expressions with the output of the COSMICS package. This allows us to 
calculate the power spectrum for the WIMP perturbations on sub-galactic scales shortly 
after matter-radiation equality (Section EJ. We consider the case of a scale-invariant 
primordial power spectrum and discuss the modihcations due to tilt and running of the 
spectral index for single-held inhationary models. Finally in Section El we estimate the 
epoch at which typical-sized inhomogeneities go non-linear for a range of benchmark 
WIMP models and examine the ehects of scale dependence of the primordial power 
spectrum on our results. To conclude we also estimate the epoch when the very hrst 
(rare) inhomogeneities go non-linear and estimate their size after collapse and speculate 
about their fate and relevance for dark matter searches. 

2. WIMPs 

Weakly interacting massive particles are generic candidates for CDM. The reason is 
that they are quite natural in extensions of the standard model of particle physics, the 
most popular example being supersymmetry. In supersymmetry models every standard 
model particle has a supersymmetric partner and in most models there is a conserved 
quantum number (R-parity), which makes the lightest supersymmetric particle stable. 
Supersymmetry models have a large number of free parameters, however in most models 
the lightest supersymmetric particle is the lightest neutralino (which is a mix of the 
supersymmetric partners of the photon, the Z and the Higgs bosons) and an excellent 
CDM candidate (see e.g. Reference (IHl)- We have focused our attention on the 
neutralino (in particular a bino-like neutralino) in our previous work P3]- Here we 
wish to dwell on the model-independent aspects of WIMP astrophysics. Our treatment 
is valid given three assumptions hold true: we assume that there is no WIMP anti- 
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WIMP asymmetry in the universe and we assume that WIMPs have been in chemical 
and thermal equilibrium with the radiation component in the early (hot) Universe (thus 
Wimpzillas are different from what we call generic WIMPs). The third assumption is 
that the elastic cross sections are dominated by exchange. This assumption is wrong 
for bino-like neutralinos, but as our previous studies have shown, the hnal result is not 
very sensitive to this detail. For simplicity, we also assume that there is only one CDM 
component, namely WIMPs. 

In the early universe WIMPs can be treated as an ideal Bose or Fermi gas 
characterised by g internal degrees of freedom, the WIMP mass m, chemical potential /i 
and temperature T. Unless the WIMPs have a net non-zero quantum number (in which 
case a mechanism is required to generate a WIMP anti-WIMP asymmetry) the chemical 
potential is negligible and, independent of the spin statistics, for T -Cm and /i -C m 
the distribution function takes the Boltzmann form and the relative relic abundance of 
WIMPs is given by: 


^wimp 2 I I 


9. 




23^/2 
TT ) 


gm 


exp (-Xcd) cosh (xcd—) • (1) 
V m / 


0 f7s|cd 

Here mpi denotes the Planck mass, H the Hubble rate, gg the number of relativistic 
degrees of freedom contributing to the total entropy of the radiation fluid and x = m/T. 
The subscripts ‘0’ and ‘cd’ refer to the present day and the epoch at which WIMP anti- 
WIMP annihilation ceased and the WIMPs chemically decoupled, respectively. In the 
following we assume /Xcd = 0. 

Equation m can be solved iteratively for Xcd: 
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After chemical decoupling, T < Ted the total number of WIMPs remains constant 
and they are kept in local thermal equilibrium by elastic scattering processes with 
relativistic particles: WIMP + L <—^ WIMP + L. As the Universe expands the 
WIMP density decreases, the elastic scattering rate decreases and the WIMPs kinetically 
decouple. The characteristic time scale between elastic scatterings is Xf = 1/Fei, where 
Pel denotes the elastic scattering rate. The average momentum exchanged per collision 
is small, of order T juiini, and the number of elastic scatterings needed to keep the 
WIMPs in local thermal equilibrium, N, is large: N ^ m/T 3> 1. The relaxation 
timescale, Xr, which characterizes the time at which the WIMP kinetically decouple, is 
given by Xj. = iVxf and is signihcantly larger than the elastic scattering timescale. 

The elastic scattering rate is given by 


Fei = ^ (vcrei)nL , (3) 

LeSM 

where is the total cross section for elastic scatterings of WIMPs and relativistic 
Standard Model fermions, ul the number density of relativistic particles of species L, 
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which are assumed to be in local thermal equilibrium and u ~ 1 in this case. The thermal 
average of cXei can be written as (agi) = where CTg* ~ {G-pm^Yw?/rri:^ sets 

the scale for the cross section and I parametrises the temperature dependence. Here, mw 
denotes the mass of the charged gauge bosons in the electroweak interaction, mz is the 
mass of the neutral gauge boson and Gp is Fermi’s coupling constant. In the Standard 
Model, elastic scattering between a light fermion and a heavy fermion is mediated by 
exchange and 1 = 0. Other channels may occur however. In supersymmetric extensions 
of the Standard Model, where the lightest neutralino is a WIMP candidate, sfermion 
exchange occurs (and if the neutralino is a gaugino, Z° exchange is suppressed), in which 
case I = 1. 

Kinetic decoupling of WIMPs happens at a temperature T^d, dehned by Tj.{TpY> = 
hf“^(Tkd). Solving this equation for Xkd = Tn/T]^d, we hnd 


^kd 


C(3) /90^7g|kd^'/' el, 


1 
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(4) 


1 /2 

where ge is the number of degrees of freedom contributing to the energy density. For 
I = 0 the kinetic decoupling temperature does not depend on the WIMP mass and is 
given by Tkd ~ 2.4 g7^^^ |kd MeV. For I = 1, Tkd ~ 34.2 g7^^^ |kd (m/lOOGeV)^/^ MeV. 

In hgure ^ we plot the dependence of the WIMP chemical and kinetic decoupling 
temperatures on the WIMP mass, for WIMPs with relic densities corresponding to the 
WMAP measurement of the GDM density (i.e. assuming that the cold dark matter is 
entirely in the form of WIMPs) lUcdm = ^cdmh^ = 0.076 — 0.156, consistent with the 
2-sigma error of WMAP P . 

While chemical decoupling happens at a temperature of around 10 GeV, kinetic 
decoupling is delayed by the large entropy of the hot Universe and takes place at around 
10 MeV. This is generic for any WIMP that is of cosmological relevance today. 


3. Cosmological perturbations of fluids 

Fluctuations in the energy density of radiation and matter come together with 
fluctuations of the scalar metric potentials. In the conformal Newtonian (or longitudinal 
or zero shear) gauge, perturbations of an isotropic and homogeneous, spatially flat metric 
are characterised by two scalar potentials 0 and Y which appear in the line element as 
(e.g. Reference j2I]) 

ds^ = a^irf) [— (1 -I- 20) (+f + (1 — 20) 5ij dxMx-^] , (5) 

where a denotes the scale factor and g conformal time. We work in this gauge since all 
sub-horizon quantities can be interpreted in terms of Newtonian physics, in particular 
the scalar potential 0 is identical to the gravitational potential in the Newtonian 
limit. This is in contrast to the synchronous gauge where the Newtonian gravitational 
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Figure 1. The dependence of the WIMP chemical and kinetic decoupling 
temperatures, Ted and T^d, on the WIMP mass (indicated by the grey scale) for WIMPs 
with Wwimp = 0.076 — 0.156. The upper and lower bands are for 1=1 (Majorana) and 
0 (Dirac) respectively. 


potential is gauged to zero. Furthermore, in the conformal Newtonian gauge there is no 
complicating residual gauge degree of freedom. 

For an isotropic fluid with energy density e, pressure P and four-velocity u^, the 
components of the energy-momentum tensor are given by 

= (e + P)u^^u, + P6^ . (6) 

Perturbations in the energy density and pressure are denoted by Se and SP respectively, 
and we treat the (small) fluid peculiar velocity Vi as a linear velocity perturbation. The 
peculiar velocity is related to the spatial component of the four-velocity by = v'^/a. 
The components of the perturbed energy-momentum tensor are given to first order in 
the perturbations by 

= -(5e , 5T; = - (e + P) v, , 5T) = 5} 5P + Pb (7) 

ct 

Here, Hb = T*- — SjT\/3 is the traceless anisotropic pressure. Quantum fluctuations 
during inflation do not seed anisotropic pressure, however once neutrinos decouple from 
the photon-lepton fluid at about T ~ 1 MeV, they build up anisotropic pressure, which 
has to be taken into account to obtain the correct normalisation of the power spectra. 
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We will include this effect below, but neglect the anisotropic stress of the neutrinos for 
all other aspects in this paper. 

In a spatially flat background it is convenient to consider Fourier modes with 
wavenumber k. For linear perturbations these modes are decoupled from each other 
and, due to the isotropy and homogeneity of the background, all modes with the same 
/c = |k| obey identical mode equations. 

We dehne dimensionless scalars to characterise the fluctuations of each fluid 

A = _ - _ 

(6 + P)a ’ 

and Ua, the modulus of the peculiar velocity, 

Ua = -*k ■ Va , 


( 8 ) 


(9) 


with the index index “a” denoting the type of fluid (e.g. radiation Pj. = or CDM 

Pcdm = 0). 

It is also useful to dehne the total energy density contrast A and peculiar velocity 
modulus v: 


a 


€ + P 


Aa 


V = 


E 


(e + P) 

e + P 


( 10 ) 


which act as sources of the energy and momentum constraints respectively. For huids, 
the pressure can be related to the energy density and the entropy by an equation of 
state Pa = Pa(ea,>S'a)- We need only consider the situation of adiabatic and isentropic 
evolution of each huid, thus 6Sa = 0, and hence 


(5 Pa = cf(5ea , 


c? = 


de^ 


( 11 ) 


“a”. 


where Ca is the adiabatic sound speed of huid 

Each of the huids individually (if they are noninteracting and dissipationless) obeys 
the general relativistic continuity and Euler equations. These equations read (here we 
have used = Pa/Ca, which is valid for linear perturbations only): 


Aa = kvs, + 3^7' , (12) 

Ua + (1 - 3Ca)7fna = “ C^/cAa - k(j) + ^/cTTa , (13) 

where the primes denote derivatives with respect to conformal time rj and P = a'/a. 
Anisotropic pressure is generated by tt which is dehned by II*/(e + P) = (5j/3 — Pkj)n. 


Any number of perfect huids can be described by one ehective imperfect huid. The 
ehective dimensionless entropy perturbation is 


5 = 


6P - ci6e 
e + P 




(14) 


where the ehective adiabatic sound speed is given by 

= w , 

^ e + P 

a 


( 15 ) 



( 16 ) 


For a Universe containing CDM and radiation fluids we And 

_ y /A A \ 2 _ ^ ^ 

“ 4 + % ^ “ 3 1 + 3|//4 ’ 

where y = a/oeq is the scale factor normalised at matter-radiation eqnality. 

We can now dehne isentropic initial conditions by demanding that 5 = 0 and 
S' = 0. The hrst condition gives Acdm = Aj. and the second, making nse of the continnity 
eqnation m, gives Vcdm = 'Vr additionally. The generalisation to an arbitrary nnmber 
of flnids is straightforward. Below we restrict onr attention to these (isentropic) initial 
conditions. 

In order to close the set of equations we also need the Einstein equations; the 
background equation 

= 47rGa2(e + P ), (17) 


and the energy and momentum constraints, which have A and v as source terms 
respectively, 

- ^ ’ ( 18 ) 

-k (V + 770) = -H')v . (19) 

Combining the spatial trace of the Einstein equations and the energy constraint one can 

write down an equation which has 5 as a source (see e.g. Reference |2I1), but we do not 

need it here. The spatial off-diagonal Einstein equation couples the anisotropic pressure 
to the metric potentials, such that 

A)2(^ - 0) = 2{n^ - n')7i. (20) 


For isotropic fluids, vr = 0 and hence ^|J = (p. In general the difference between ip and 0 
can be neglected on subhorizon scales {k^ P.). 

To normalise the spectra of cosmological perturbations, we have to connect solutions 
in the radiation, matter and dark energy dominated epochs with the primordial 
fluctuations that are produced during the epoch of cosmological inflation. To do that 
it is most convenient to introduce a quantity that is independent of the choice of the 
hypersurface of constant time and is constant for modes in the superhorizon regime 
{k <^P). Following Bardeen j22| we dehne 

c = A/3-0, ( 21 ) 


which is the curvature perturbation on uniform density hypersurfaces or the density 
contrast on uniform curvature hypersurfaces. We see immediately from the continuity 
equations m, that this quantity could be dehned for each huid or for one effective huid 
and that it is constant for regular solutions on superhorizon scales (assuming that there 
are no entropy perturbations) [22] • In the following we normalise the power spectrum 
of C to the value measured by WMAP, where this variable is denoted —TZ. 

During radiation domination (cr S> ecdm), the perturbation equations ([T^ and m 
can be solved exactly for the cold dark matter and isotropic radiation (vr = 0 ^ 0 = 0) 
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fluids, on all scales 1201. In the conformal Newtonian gauge we find for the dominant 
growing mode 


[K 


= -2Co 


JiW 


K 


. rad . X „ ^ 

A,, {k) = 3Co 


Jo(k) 


2-h kji\k) 

K 


= \/3Co [Kjo(K) - 2ji(K)] , 

^rdm(«) = 6Co 

1 - io(« 


ln(K) + jo{K) - - Ci(K) + 7 e - ^ 

K 2 




= 2^3(0 


K 


( 22 ) 


(23) 


where n = krj/\/3, (^o is the value of ( in the superhorizon limit, 7 e Euler’s constant, Ci 
the cosine integral and jn are the spherical Bessel functions, all dehned as in 123 . It is 
easy to check that this is the isentropic mode, since in the superhorizon limit Acdm Aj. 
and Ucdm ^ Vr. 

On superhorizon scales {k -C H), 0 ~ —2(o/3, ~ A — ik/H)Co/3. On 

subhorizon scales (fc S> Td), a hrst order expansion of the exact solutions, (1^ and 
in K gives 

cos(k) 


\K) 


-2Co 




~ — 3CoCos(fi:), ~ ■\/3(^osin(K) , 
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^cdm(«) - 6C0 


ln(K) + 7 e - 


rad , . 2^/3Co 

^'cdml^) - - 

K 


(24) 

(25) 


During radiation domination, the radiation density and velocity perturbations on 
subhorizon scales oscillate with constant amplitude, thereby generating a Newtonian 
gravitational potential that decays like 1 /k^, while the CDM density perturbations 
grow logarithmically and the velocity perturbations decay with 1/k. From (12411 it can 
be seen, that the Newtonian gravitational potential acts like an external held on the 
evolution of CDM density perturbations in this regime. 

In the next section we discuss the viscous coupling of WIMPs to radiation and its 
effect on the evolution of WIMP density perturbations during radiation domination. 
While the viscous coupling has important consequences for the spectrum of WIMP 
perturbations, the evolution of the radiation huid is unaffected. 


4. Kinetic decoupling 

Elastic scattering processes around Tkd lead to viscous coupling of the WIMP and 
radiation huids, which results in collisional damping of WIMP perturbations, as shown 
by two of us in |I2|. The WIMP perturbations disperse due to bulk and shear viscosity. 
The strength of these damping mechanisms is described in terms of local transport 
coefficients Cvis and r^vis for bulk and shear viscosity. Close to local thermal equilibrium. 
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Cvis ~ S/SnTrreiax and rj^is ~ nTr^eiax. The WIMP sound speed in that regime is given 
by -y/^hTsy^TT?^. The WIMP perturbations that survive the viscous coupling prior to 
kinetic decoupling are then the initial conditions for the free streaming regime. 

We now calculate the spectrum of WIMP density perturbations on the surface 
of kinetic decoupling by solving the linearised Navier-Stokes equations. For a more 
general and complete derivation, see ^21- Dissipation is a subhorizon phenomenon. The 
linearised Navier-Stokes equation for the WIMP fluid on subhorizon scales {k S> 'H), is 
given by 


V • + 

^wimp ' 


Cvis + 4/3r7^is 


’^^wimp *^wimp^^wimp5 


-Wimp 


(26) 


where a gravitational forcing term and energy dissipation due to the large-scale 
inhomogeneities of the radiation fluid, as well as horizon scale contributions, have been 
dropped. Similarly, the continuity equation m becomes = /cUwimp- These 

equations can be combined to give a second order differential equation for the evolution 
of 


Cvis + 4/3?7vis k^ 2 ^2 ^ . -0 ('27') 

^wimp ' ^wimp ' ^wimp ^ ^wimp ' J 

^wimp ^ 

Note that the final term should not be neglected as its coefficient (cPmp ^ T/m) is of 
comparable magnitude to that of the dissipative term [(T/m)(rreiax/^)] for t ~ Treiax- 
This is the equation of motion of an oscillator with time dependent friction, due to 
the bulk and shear viscosity, and hence A^imp oscillates with complex frequency. The 
real part of the frequency is proportional to the isentropic sound speed in the WIMP 
fluid and describes the propagation of perturbations in this fluid, while the imaginary 
part describes the damping of the perturbations caused by the transfer of energy and 
momentum, due to the viscous coupling, from the WIMP fluid to the radiation fluid 
(which acts as a heat bath). 

Let US stress that our treatment is restricted to A; 3> 7Y. Near the horizon there are 
extra forcing terms on the right hand side of (EH). There are two types of such terms. 
Firstly, there is the gravitational pull of the oscillating radiation fluid. Secondly, as long 
as Cwimp ~ Crad the oscillatious of the radiation fluid can drag the CDM fluid along by 
means of multiple scatterings in a preferred direction. As soon as the relaxation time 
scale exceeds the oscillation period of a given mode, this mechanism can no longer be 
active. This coherent effect is not accounted for by our viscosity terms, which reflect 
the dissipation due to the homogeneous component of the radiation fluid. There are 
also some extra terms from the Hubble expansion on the left hand side of (EZD- 

The numerical treatment of Loeb and Zaldarriaga ^H] includes these terms. Their 
results conhrms that the essential features of the subhorizon damping are captured by 
our calculation. They also show that the additional terms play an important role in 
determining the precise position of the maximum of the CDM power spectrum. We 
agree with their conclusion that an accurate calculation (better than 10% accuracy) 
must be based on a numerical computation including all of the terms. 
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The WKB solution to (j27|l for the envelope of the WIMP density oscillations is: 

Cvis + 4/3?7vis 


^wimp(^; ^kd) ^wimp {k, r]i) exp 


'Vi 


2 c wimp 


-dr; 


(28) 


where Awimp(fc, vO is the initial primordial density perturbation and r;kd is the conformal 
time at kinetic decoupling. The damping term can be written as 


Da{k) = 


_ ^wimp(^)hkd) 


^wimp(^) hi) 


= exp 


(29) 


where fed is given by 
kd = 


Cvis + 4/3r;vis^^^ 
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(30) 


This scale corresponds to a length scale ~ 10“^/if at kinetic decoupling. The total 
WIMP mass contained in a sphere with radius rr/fcd is Md ~ IO^^^Mq. Here we 
assumed a WIMP mass typical for the lightest supersymmetric particle and a typical 
kinetic decoupling temperature for a weakly interacting particle of such a mass nacii. 

Once the WIMPs have kinetically decoupled, they can be described as a separate 
fluid. On scales larger than the free streaming length, we can approximate it as a 
pressureless fluid. Since kinetic decoupling happens well before the end of radiation 
domination the gravitational evolution of the WIMP density perturbations on arbitrary 
scales during radiation domination -C er) is given by (OHl) . 


5. Free streaming 


After kinetic decoupling the evolution of the WIMPs at the smallest scales is described 
by the distribution function /(x, gn, r;), where q and n are the norm and direction of 
the comoving 3-momentum, i.e. gph = q/a. The distribution function is governed by the 
collisonless Boltzmann equation, which reads in a flat Friedmann model (see e.g. I2SI) 


f d dx* d dq d dn,- d\ 


0 . 


(31) 


In local thermal equilibrium we have /(x, gn, ?;) = /o(g, ?;), which in the case of 
non-relativistic particles becomes the Maxwell-Boltzmann distribution function (with 
g internal degrees of freedom). 


fuM = jijs exp (^ , exp 


{q/ckf l‘2.m 
T~ 

Wimp 


(32) 
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Here, T^imp denotes the temperature and /i the chemical potential of the WIMPs. Both 
depend on time, but /mb(?) does not. 

In the next step we have to investigate deviations of the distribution function of the 
WIMPs in order to describe an inhomogeneous universe. We will do so by considering 
small perturbations in the Boltzmann equation away from the local thermal equilibrium 
solution. Close to local thermal equilibrium we can write 

/(x, gn, T]) = /o(g, t]) + ^/(x, gn, t]) . (33) 

The last term in m does not contribute at hrst order in the perturbed quantities. 
This is because drij/di] and df /duj are individually hrst order terms; in the absence of 
metric perturbations free-streaming particles do not change direction and the zero-th 
order distribution function /o(g, "g) is independent of direction n. Let us now turn to 
the third term in (jSH). From the geodesic equation it follows that dq/dr] is proportional 
to derivatives of metric perturbations. Free streaming only occurs on scales well below 
the horizon (due to the small velocities of WIMPs) and on subhorizon scales the metric 
perturbations (j) and '0 are negligible — consequently the third term in m can also 
be neglected. The collisionless Boltzmann equation for subhorizon scales can then be 
rewritten in terms of the spatial Fourier transform of 5/(x, gn,-g), which we denote by 
5/(k, gn,g): 

+ ^ — = (34) 

\ori m J 

Note that the second term in this equation depends on the direction n of the comoving 
momentum only via its angle with respect to the wavenumber k. 

Before solving equation (El the initial perturbations have to be specihed. At 
kinetic decoupling we can assume that T = T^imp and we drop the index “wimp” in 
the following. The WIMP phase space distribution around kinetic decoupling is close 
to a Maxwell-Boltzmann distribution El- The WIMP density perturbations (with 
7^ -C A; -C fed) which are present on the surface of kinetic decoupling need to be taken 
into account. To do so we consider small fluctuations ST and d/i of the thermodynamic 
variables T and n in and expand the resulting distribution up to hrst order: 

f (^.-M) = ) + (? + T + ® (^7 . (35) 

As we are interested only in effects of hrst order in the perturbed quantities, the adiabatic 
equations of motion can be used to establish the relationships between temperature and 
chemical potential huctuations as well as WIMP perturbations on the surface of kinetic 
decoupling. From the conservation of entropy per WIMP we hnd 6n/n ~ {3/2)ST/T 
and thus ST/T ^ 2Awimp(A(,'gkd)/3. Inserting this in the Gibbs-Duhem relation 
S^n/T) = SP/{nT) — (e -I- P)/{nT){ST/T) and using the previous relation, we hnd 
S{fi/T) ^ — (2m/3T)Awimp(A:, gkci), so that the deviation from local thermal equilibrium 
at kinetic decoupling is given by 

(T\ (k,q,ni) = + 0(g‘) . 


(36) 
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(/c,g,?7kd) exp [-il{q,T]) k 


n 


(37) 


This has an obvious interpretation, the fractional perturbation of the distribution 
function is proportional to the ratio of the kinetic energy of an individual particle with 
comoving velocity q/m and the thermal energy 3T/2 times the density contrast. 

The solution to the free streaming equation ()34|1 is then 

I) (I 

where I = liq,^]) is the comoving distance a WIMP can travel freely in the background 
space-time during the time interval rj — ?7kd: 

We can evaluate that this integral for the matter-radiation universe, in which we can 
write the scale factor as 

2 


a{r]) = a. 


eq 




+ 




Using the previously introduced notation y = a/oeq some algebra dually yields. 


= — 2 /kdln 
m 


l/kd 


2/ + 2 (1 - vr+^) 




n 


(39) 


(40) 


eq 


y 2/kd -|- 2 (l — \/l -|- l/kd) 

where gkd = <?/akd is the modulus of the physical momentum at kinetic decoupling. 
As 1 / -C 1, the physical distance of travel of free streaming particles approaches 
[(<?kd/''^)l/kd]'\/2ln(4/j/kd)[a/7feq], the hrst square bracket is the physical velocity at 
kinetic decoupling and the scale is set by the physical size of the Hubble horizon at 
equality (the other term in square brackets). 

The free streaming of the WIMPs generates a further (collisionless) damping 
mechanism for WIMP density perturbations. After averaging over the Maxwellian 
distribution of velocities, we can estimate the damping scale by replacing the physical 
velocity of an individual WIMP by the mean velocity (thermal velocity) of WIMPs at 
kinetic decoupling. Thus the comoving length scale of free streaming 

with hkd = \/^gT/m. In order to characterise the spectrum of perturbations which 
survive collisionless damping we introduce the comoving free streaming scale fcfs oc l//fs. 
More precisely we dehne 

v/2 


kfsiv) = 


(42) 


^rkd/m[/(gkd,h)/(<?kdM)] ’ 
where the numerical prefactor is chosen for later convenience. With the help of ii we 
hnd that the comoving kfg becomes approximately constant as a/oeq ^ 1 and is given 
by 


kf: 


f) 

-Lkd/ 


1/2 


®eq/®kd 


eq 


ln(40eq/®kd) 




eq ) 
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1.70 X 10® (m/100 GeV)i/2(Tkd/30 MeV)!/^ 


(43) 


Mpc 1 + ln(rkd/30 MeV)/19.2 

This expressions depends on only throngh the logarithm, we therefore set 

it equal to the WMAP best fit value, (Um = 0.14 [Tj. The corresponding length scale at 
matter-radiation equality is ~ 10“®/i7 and the total WIMP mass contained in a sphere 
of radius 7r//cfs is Mfg ~ IO^^Mq. 

The WIMP density contrast is related to the distribution function by: 

/dG 5/(k,gn,r/) 


A 


Wimp 


(k, v) = 


(44) 


47r/“dgg2/o(g,??) 

with dqq'^dfl denoting the volume measure in spherical coordinates in q-space. The 
denominator in (I44j) is the comoving mass density of CDM and the numerator is the 
CDM mass fluctuation due to the deviation from local thermal equilibrium. Using (Eni) 
and dSZl) in (I44|l . we find after integration over 

1/2 


A. 


Wimp 


(k, v) = 


TT 


A 


- 3/2 


X exp 


wimp (k, 7/kd) {mT]^d) 

sin iliqi,d,v)k) 


dqkdq^ 




kd 


3mTkd 


Qkd 


2mT 


Kqkd,ri)k 


(45) 


The remaining integration can be found as expression 3.952(5) in Reference 

As a result of the integration we find the suppression of the WIMP density contrast 
due to free streaming to be 


D{s{k,v) = 


A^ijjip (k, 7/) 


Awimp(k, 7/kd) 


1 - 


A 

kfs 


exp 


A 

fcfs 


(46) 


We see that the free streaming of WIMPs results in exponential damping of the 
WIMP density contrast, similar to that produced by collisional damping (justifying 
the terminology ‘collisionless damping’ for the effects of free streaming). There is one 
difference, however; the ratio of the WIMP kinetic energy to the thermal averaged 
kinetic energy in dsni) leads to a polynomial pre-factor. This expression is valid for 
k/kfs < 1, as terms of order {k/kisY have been neglected in the polynomial. 

The net damping is the product of the collisional (due to viscous coupling to the 
radiation fluid) and collisionless (due to free streaming) damping terms [(jSHl) and (HHll 
respectively], D{k) = Dd{k)Dis{k)\ 


D{k) 


Awimp(^) ^) 
Awimp (^) ^i) 


1 


2 

3 



exp 



. (47) 


Since fcfs <C kd, the cut-off in the power spectrum is determined by the free streaming 
scale kfs- 

In figure |2 we plot the variation of the characteristic damping and free streaming 
comoving wave numbers with WIMP mass, for WIMPs with relic densities corresponding 
to the WMAP measurement of the CDM density (i.e. assuming that the cold dark 
matter is entirely in the form of WIMPs) oicdm = = 0.076 — 0.156, consistent 
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Figure 2. The variation of the characteristic damping and free streaming comoving 
wavenumbers, fcj and kfs, with WIMP mass (indicated by grey scale). The lower and 
upper bands are for 1=1 (Majorana) and 0 (Dirac) respectively. 


with the 2-sigma error of WMAP P . For concrete calculations we will concentrate on 
four benchmark models which span the range of most plausible WIMP properties. The 
details of these benchmark models, including the values of fed and kfs, are tabulated in 
Table [T] Models B and C are very close to the bino-like neutralino cases of our previous 
work ^3]. Models A and D show that there is more spread of the predicted damping 
scale if the strong assumption of a bino-like WIMP is dropped. 

In order to compare our analytic results with the numerical result of reference jl8j . 
we express the cut-off scales in terms of matter mass enclosed in a sphere of radius 
R = vr/fcfs, Mcnt{R) ~ 4.9 x 10“®MQ(l/fcpc)^. For a decoupling temperature of 25 
MeV, we obtain from pH] a mass of 6.4 x which is to be contrasted with 

1.5 X 10~^Mq from our calculation. So, the difference between our analytic result and 
the numerical result including effects near the horizon scale at decoupling is a factor of 
a few, instead of two orders of magnitude as claimed in US]. 

The approximate equality of the kinetic decoupling and free streaming scales (to 
within an order of magnitude) can be understood by comparing the corresponding 
physical length scales at equality (after equality the comoving free streaming scale only 
grows logrithmically). The free streaming length ()41|1 is given by the product of the 
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Table 1. Benchmark WIMP models. 


Ref. 

/ 

m (GeV) 

Ted (GeV) 

Tkd (Mev) 

fed (pc"^) 

fcfe (pc“^) 

A 

0 

100 

3.6 

1.6 

14 

0.42 

B 

1 

50 

1.9 

21 

39 

0.94 

C 

1 

100 

3.7 

25 

61 

1.5 

D 

1 

500 

17 

37 

180 

4.0 


WIMP velocity at equality, the logarithmic growth factor and the Hubble radius at 
equality: 

(48) 

^eq ^kd 

From pi|) . the collisional damping length at kinetic decoupling is given by the square 
root of the product of the viscosity terms with the time of kinetic decoupling, divided 
by the WIMP mass density. The viscosity terms are proportional to the product of the 
WIMP mass density, the WIMP sound speed Cwimp squared and the relaxation time. 
At kinetic decoupling the relaxation time is by definition of order the physical Hubble 
radius. The collisional damping length at equality is thus 

~ Cwimp.Rkd ~ Cwimp Req ; (^9) 

(2kd ^eq 

where we have used H cx 1/a? during radiation domination and Cwimp ~ ■\/Tkd/m. At 
kinetic decoupling the sound speed and the average WIMP velocity are approximately 
equal and therefore /fs ~ Zd In (oeq/akd), i-e. the length scales are roughly equal at 
equality. 

This can also be seen from the Reynolds number of the CDM fluid which is 
given by Re = 2RHC^impP/v rsj {Rh/T) c^i^^m/T. At kinetic decoupling we have Re 
^/mjT 1/cwimp- Thus the sound speed of the CDM (equivalent to the mean 
particle velocity in the non-relativistic case) at kinetic decoupling plays a fundamental 
role in both processes. Kinetic effects dominate over friction once Re is large. Thus for 
Tkd ~ 10 MeV and m = (9(100GeV) we hnd Re ~ 100. We note that by fluid engineering 
standards, this is still a pretty small Reynolds number. From this consideration we can 
also see that kinetic decoupling is a very fast process, as Re ~ (for the case I = 1 
in Eq. Ej). 

6. Evolution of CDM perturbations 

Now we turn to the gravitational growth of CDM perturbations. So far we have 
calculated the damping of the WIMP density perturbations due to collisional and 
collisionless damping (HID, and the gravitational evolution of CDM perturbations for 
^ ^cdm; (El- What remains is to hnd a solution for the gravitational evolution of 
CDM perturbations during the matter and dark energy dominated epochs. In the 
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following we find analytic approximations to the radiation-matter and matter-dark 
energy transitions, which hnally can be smoothly joined with (I2d|l to provide the correct 
normalisation to CMB measurements. 

Since we are interested in the subgalactic scales only, we can restrict our attention 
to the sub horizon modes when Ccdm becomes comparable to e^. We include neutrinos in 
the radiation component in order to allow an analytic treatment, i.e. their anisotropic 
stress is neglected. This leads to errors of around 10 % We also neglect the baryon 
inhomogeneities. At early times the baryons are tightly coupled to the radiation fluid, 
and photon diffusion damping rapidly erases small-scale perturbations in the baryon 
fluid at z ^ 10® to 10®. On small scales the tight coupling breaks down prior to 
recombination, and the baryon perturbations grow, however Ab Acdm still OH] . 
Post decoupling on scales fc > fcb ~ lO^Mpc”^ the residual electrons allow transfer of 
energy between the photon and baryon fluids so that thermal pressure prevents the 
baryon perturbations from growing, until z\, ~ 150 EH) • As we are interested in 
CDM perturbations on small scales at early times, we can neglect the perturbations in 
the baryon fluid. 

For scales which enter the horizon sufficiently long before matter-radiation equality 
the logarithmic growth of CDM perturbations provides a situation in which, after some 
time, CmAm ^ e^Ar (i.e. only the CDM perturbations are important) even during 
radiation domination EH- We also have CmAm 3> CdeAde, where the index de stands for 
dark energy. In the following we will assume a cosmological constant, for which Ade = 0. 
Thus we only need to keep Am in A on the rhs in the Poisson equation (HED- 

With these two approximations (subhorizon scales and cdm fluctuations dominating 
as the source in the Poisson equation) we can simplify the equations (fT^ and (fT^ to 

^cdm ^Ucdm; ^cdm T "^Ucdm 5 (^0) 

and the Poisson equation reads 

- P(j) = dvrCa^ecdmAcdm • (51) 


It is now most convenient to combine these equations to a single one and to use the 
scale factor instead of conformal time: 


.d^A, 


cdm 


da^ 


+ 


1 - 


P\ dA 


cdm 


3 ^cdm 


da 


A, 


cdm 


= 0 . 


(52) 


6.1. Radiation-matter transition 

For the radiation-matter transition (at which point the dark energy, or cosmological 
constant, is negligible) simplihes further 

1/(1 + y)-^^cdm + + \y^ ~ ~ /b)Acdm = 0 , (53) 

where y = a/oeq = Cm/Cr and /b = cab/wm is the baryon fraction, with best £t value from 
WMAP ff, = 0.17 PJ. The exact solution to this equation is a combination of Legendre 
functions of the hrst and second kind: 


Acdm(fc, y) — + y) + i? 2 (^)C..(\/l + y). 


(54) 
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with index v^fb) = ( 1/25 — 24/^ — l)/2. 

We now determine the constants Bi 2 {k) by matching the y <C 1 expansion of 
equation dSl to the subhorizon limit of the radiation domination solution (jsni)- We hnd 


Bi{k) = 6 C 0 


In 




+ h 


B2{k) — — 12(^0 


where fceq = "^eq and 


Kh) = 2 In 


— 7e 


1 

2 


V 


2r>] 

rM 


(55) 

(56) 

(57) 


where T'[v\ is the derivative of r(z/) with respect to v. 

Finally expanding (ED for y ^ 1 we hnd that during matter domination, for scales 
that enter the horizon signihcantly before matter-radiation equality and z < Zh 


^cdm{y) = 6Coc(z/) 


In I A 


eq 


+ h 


where 


c z/ = 


r[i -|- 2z/] 


(58) 


(59) 


2^r2[i + i/] ’ 

e.g. z/(/b) = 1.79 (2), c[z/(/b)] = 1.36(3/2) and 5(/b) = -1.56 (-1.74) for /b = 0.17 (0). 
Note that before ;2b, CDM density perturbations grow as Acdm oc Later the baryons 
follow the CDM and the matter huctuations grow as a. For the peculiar velocity and 
the Newtonian gravitational potential we obtain 


'^^cdm(2/) 

0 ( 1 /) 


'^eq 

k 



Acdm(2/) 




3 / fceq A 
A\k ) 


(1 - A) Add 

y 


(60) 


In the following, we omit the subscript ‘cdm’. 

For redshifts z^q^ > z > z-^, (between matter-radiation equality and the epoch at 
which small-scale baryon perturbations start growing) we hnd the transfer function for 
the CDM density perturbations for modes which satisfy fc > fcb 

n 2 


TA{k,z) = {Qcf 


In 


k 


k 


+ b 


eq 


1 + 7 


eq 


1 + Z 


(61) 


and the transfer function for the Newtonian gravitational potential on these scales is 
given by 


T^{k,z) = 


27(1 - /b)c 


In 


k 


+ b 




1 + 7 


eq 


1 + ^ 


u-2 


(62) 


The transfer function for the velocity depends on the initial time and is therefore not a 
very useful quantity. 
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6.2. Normalisation including anisotropic stress from neutrinos on superhorizon scales 


Above we have made the assumption that (f = i/j as we neglect the anisotropic stress 
of neutrinos, i.e. tTj, = 0. For the normalisation this leads to unacceptably large errors. 
We therefore take the superhorizon anisotropic stress of neutrinos into account in the 
normalisation. 

Let us dehne the radiation fraction in neutrinos 

a = — = 0.405, (63) 

after e’''e“ annihilation and for three massless neutrinos. The regular solution for 
neutrinos at superhorizon scales during the radiation dominant epoch gives isniEi 





^rO — 


A 


cdmO 



(64) 


and thus 


Co 



(65) 


The conservation of ( on superhorizon scales then implies that, for z < ^eq but before the 
start of dark energy domination, the normalisation of CDM fluctuations on superhorizon 
scales is given as (now we can again assume <p = since the neutrino anisotropic stress 
is now suppressed by 

9 / 4q;\ 

Acdm = --U + —j0o. (66) 


When normalising to Co, which is provided by the WMAP measurement (encoded as 
the variable A), we do not have to take special care as the WMAP value includes the 
effect of neutrinos on superhorizon scales. However, when comparing our analytic result 
with the results of the COSMICS code below, we have to include this correction, as the 
initial condition in the code is (/)o = — 1. 


6.3. Accuracy of approximations 

In Figure El we plot the CDM density contrast at = 300 and 100 as a function 
of comoving wavenumber for the WMAP best £t total matter and baryon densities, 
<^cdm = 0.116, (Ub = 0.024 (/b = 0.171 and 6(0.171) = —1.562), and for zero baryon 
density, oJcdm = 0.14, cUb = 0.00, with h = 0.72, using our analytic expressions and also 
using the “lingercon” Boltzmann solver from the COSMICS package (which includes 
massless neutrinos) |SSI- To allow direct comparison with the output of COSMICS we 
take initial conditions (during radiation domination) here such that Acdm = 3/2 (i.e. 
00 = — !)■ The rapid oscillations of Ar on sub-horizon scales means that it is not 
possible to accurately numerically evolve the coupled perturbation equations for the 
sub-galactic scales {k 3> lO^fceq) that we are interested in. We see that our analytic 
expression accurately reproduces the shape of A for k > lO^fceq and the normalisation 
is accurate to ~ 10%. The accuracy initially increases as 2 ; becomes smaller, which can 
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log,o(k/kj 


Figure 3. The CDM density contrast A at 2 ; = 300 (lower lines) and 2 : = 100 (upper 
lines) for Wcdm = 0.116, Wb = 0.024 and Wcdm = 0.14, Wb = 0.00 (bottom and top 
line of each pair respectively), using our analytic expression itKHll (solid line) and also 
the COSMICS package (dotted line). The analytic superhorizon normalisation I 66 II is 
shown as a dashed line for k < lO^^fceq- The vertical dashed line denotes kb- Our 
approximation applies best at fc > fcb but, as the comparison shows, it is also pretty 
good for smaller wavenumbers. 


be easily understood from the fact that we assumed that 2; Zeq- Our calculation then 
becomes less accurate for z < Zh- 

We also see that baryons have a significant effect on the growth of the density 
contrast; a.t z = 300 A is roughly 40% smaller in a Universe with the WMAP best £t 
energy densities (oJcdm = 0.116, cUb = 0.024) than in a Universe with the same total 
matter density, but no baryons ((Ucdm = 0.14, cUb = 0.00). 

We conclude that our analytical accuracy is good enough in the light of the present 
uncertainties in the cosmological parameters and primordial power spectrum. 

6.4- Matter-dark energy transition 

Let us stress before going into any details here, that for the subgalactic scales we are 
interested in, the following calculation is not directly relevant, since those modes become 
non-linear long before the onset of dark energy domination. However, we need to take 
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the suppression of perturbation growth into account on larger scales when we make 
contact with the measurement of ag below (section |H1). 

In ACDM cosmologies at late times the cosmological constant, or dark energy, 
dominates the Universe, leading to accelerated expansion and the suppression of the 
growth of density perturbations. For this situation equation (IS2D can be simplihed 
neglecting the presence of radiation. We dehne the scale factor relative to that at the 
epoch at which the matter and cosmological constant densities are equal: u = a/a^cfl 
and aeq 2 = oo and thus u{ao) = The evolution of the matter 

density contrast is now governed by 

l + 2n3^ d . 3/1 


u 


dn2 “ 


3 

+ 2 


d 

U — lXr 
au 


1 + J du 2 \1 + 

which can be obtained from (IS2D using the change of variables Am = uw and = —z. 
The function w{z) obeys a degenerate hypergeometric differential equation. Its general 
solution is [see jHIj, sec. 2.2.2 case 2, eqs. 2.9(1) and 2.9(18) with a = 1,6 = 1/3, c = 
11 / 6 ] 


A„ = Ci(k)u2F,(l, i; A; -u^) + 1 + 




( 68 ) 


The hrst solution is regular for small u, the second one is singular. We have to 
join the regular solution to the growing mode (EHD. In the limit u —0 we have 
2 ^ 1 ( 1 ,1/3; 11/6;—M^) ^ 1 and thus Ci{k) = (3/2)i?i(/c)(l + ;2eq)(f^A/f^m)^^^ and 
C'2(fc) = 0. 

We can now put everything together to obtain the linear matter density contrast 
today 


^m(fc) — 9Co(l + ^eq) 


In 


A 

k 


eq 


+ b 




1 . n. 

3 ’ 6 ’ 


Ua' 
ilm ' 


(69) 


An excellent analytic approximation to evaluate the hypergeometric function for 
f^A/f^m > 1 is given by the asymptotic expansion (see j21I) 


M2-Fi(i,4;if; -«“) 


2r(|)r(A) 


TT 


1 \ 5 1 , 


^ 1.25 0.719 1 , 

1.437- — H- - -h C>(—) 

u° 


(70) 


which has an accuracy better than 1% for Um < 0.4. 

The growth of structure comes to an end asymptotically; in a pure CDM model 
the density contrast grows as u, however in a ACDM universe it only grows by a factor 
1.437 after matter-dark energy equality. An equivalent solution has been obtained by 
Eisenstein jSB] in terms of elliptic functions, however he expanded the elliptical functions 
for small u, which would be more appropriate for flA/Um < 1. One could also write 
the above hypergeometric function (or elliptical functions of Eisenstein) in terms of 
associated Legendre functions, but it seems to us that nothing is gained by doing so. 
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7. Linear power spectrum 

In this section we present the linear dimensionless power spectra (defined as Vx{k, z) = 
{k^/2n'^){\X{k, z)\'^)) normalized to the WMAP measurements ffll IHH]!. 

7.1. Scale invariant spectrum 

For simplicity we start with a scale-invariant primordial power spectrum and assume 
that gravitational waves have a negligible contribution to the CMB anisotropies. We 
find for k > kh and 2 ;eq S> 2 ; ^ z\, 



where A = 0.9 ±0.1 at the pivot scale ko = 0.05/Mpc according to Reference [Tj. 


Note that from the WMAP data the spectral index n = 0.99 ± 0.04 is consistent 
with the scale-invariant Harrison-Zel’dovich spectrum n = 1. The scale of equality 
is fceq = (0.01/Mpc)(a;m/0.14) and 1 ± Zeq = 3371(a;m/0.14). 

Figure 0] shows the power spectrum for the WIMP density contrast at a redshift 
of 300, close to the end of the linear regime of structure formation, with and without 
the effects of collisional damping and free-streaming for the four benchmark WIMP 
models introduced in SectiorO In accordance with the WMAP best fit values we 
take cUcdm = 0.116 and /b = 0.17, and we assume a scale invariant primordial power 
spectrum (i.e. n = 1). It can be observed that the cut-off of the power spectrum is 
indeed very sharp (with a maximum close to the cut-off). As discussed in section |S] the 
characteristic free-streaming wave-number is smallest for model A (Dirac WIMP with 
1=0) and increases with increasing mass for the Majorana WIMPs (models B-D). This 
is clearly reflected in the position of the cut-off in the power spectrum. 

7.2. Scale dependent spectrum 

On the scales probed by the CMB [0(0.01 — 0.1)Mpc~^] the primordial power spectrum 
is close to scale invariant PP. The free-streaming scale fcfs ~ 10®Mpc“^ is seven orders of 
magnitude smaller and hence even a very small scale dependence of the power spectrum 
could significantly change the amplitude of the power spectrum close to the cut-off, and 
hence the red-shift at which the first WIMP halos form. 

To assess the effects of possible scale dependence of the primordial power spectrum 
we consider three benchmark infiation models which span the range of possible power 
spectra for simple infiation models: a. V = chaotic infiation model, power law 
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logio(k Mpc) 


Figure 4. The dimensionless power spectrum of the WIMP density contrast at z = 300 
for our four benchmark WIMP models assuming a scale-invariant primordial power 
spectrum (full lines, from left to right models A, B, C and D). Without the effects 
of collisional damping and free streaming, the power spectra would be given by the 
dotted line. The vertical dashed line denotes kb, the wavenumber below which baryons 
follow CDM. Our approximations are optimised for k > k\,. 


inflation and a hybrid inflation model. These models span an interesting region of 
inflationary parameter space, for a more detailed discussion see ffl. 

There are also uncertainties in the parameterization of the power spectrum. The 
most commonly used parameterization is 


V{k) = V{ko) 


n{ko)-l+^a(ko)\n{k/ko) 


(74) 


where a{k) = dn/dk. An arguably more appropriate parameterization over a wide-range 
of scales is gni^ 


V{k) , fk 

—— — oo -b ai In I — 
V{ko) \ko 




(75) 


A small difference between the two parameterizations can be used as an indicator, that 
the slow-roll approximation is justihed for the model at hand m- 

The spectral index, n, its running a and the alternative expansion co-efficients 
ttn depend on the inflationary potential and are most conveniently expressed in terms 
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of the horizon flow parameters EH which are dehned as: eo = H{Ni)/H{N), where 
N = ln(a/aj) is the number of e-foldings of inflation since some initial time ti, and 


_ din |e„ 
^”+1 = dAT 


n > 1 . 


(76) 


The horizon flow parameters are related to the traditional slow roll parameters, 
e = (mpi/167r)(l/7^)^) V = ("^pi/ 8 ^)(^"/^) = (mpj/87r)^(l/V”7^^) where 

' = d/d 0 , as ei = e, €2 = 2 e — 2ri and 6263 = 4e^ — 6er] + 2 ^^. 

To hrst order n — 1 = —2ei — 62 (see e.g. Ref. [32] for the higher order terms) and 
a = —26162 — 6263 IH]. oq, oi and 02 are given by equations (26) - (28) of Ref. jl2]- 


7.2.1. chaotic inflation In this model 61 = 62/2 = 63/2 = l/( 2 At-|- 1 ) (see e.g. 

Ref. ^ 3 ]) and we take N = 55 [321 • This gives, including the higher order terms in the 
expression for n, n — 1 = —0.03611 and a = —6.49 x lO”'^. 


7.2.2. Power law inflation In power law inflation the scale factor grows as a oc 
(with p > 1) and 61 = 1/p with all other horizon flow parameters equal to zero. We pick 
p = 55.4 so that the spectral index is the same as for the nPcjfl chaotic inflation model 
(n — 1 = —0.03611). In this case there is no running of the power spectrum (ck = 0), 
i.e. the spectral index is constant. 


7.2.3. Hybrid inflation Our benchmark hybrid inflation model has potential 

2 " 


V = V. 




mpi 


(77) 


and we assume that the hrst, false vacuum term, in the potential dominates so that 
"C 62 , with 62 constant. We take 62 = —0.014, the 2-cr lower limit from WMAP and 
2dF found in Ref. ^3] giving the largest increase in the density contrast with increasing 
fc, so that n — 1 = 0.03550 and a = 0. 

The primordial power spectra of these three models are plotted in Figure El for 
both parameterizations of the power spectra (equations (ITDl and (Ei). We see that 
the amplitude of the primordial power spectra on the free-streaming scale lO^Mpc”^ 
varies by a factor of ~ 2.5 (equivalently the amplitude of the huctuations varies by 
~ 7^ 1.6). We also observe that the two parameterisations do not give rise to a 

signihcant difference, thus we will stick to the more traditional power-law shape in the 
following. 

In hgureiniwe plot the processed power spectra for these primordial power spectra 
and also a scale invariant primordial power spectrum for WIMP benchmark C, with 
cosmological parameters hxed to the WMAP best £t values as before. The variation 
in the amplitude of the primordial power spectra for k ~ 10 ®Mpc“^ translates directly 
into a variation in the peak amplitude of the processed power spectra. 
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Figure 5. The primordial power spectra of the three benchmark inflation models 
discussed in the text (from top to bottom: hybrid, power law and chaotic 

inflation) for the standard power law parameterisation of the power spectrum (equation 
solid line) and for the expansion in ln(fc/A:o) (equation (I75II . dotted line). 


8. The first structures 

The collisional damping and free streaming of WIMPs lead to a cut-off in the WIMP 
power spectrum, which sets the typical scale for the hrst halos in the hierarchical 
picture of structure formation. We estimate the redshift at which typical fluctuations 
on comoving scale R go nonlinear via 

cr(i?, Zni) = 1, (78) 

where cT(i?, z) is the mass variance dehned by 

dh 

a\R,z)= / W\kR)V^{k,z)—, (79) 

where W{kR) is the Fourier transform of the window function divided by its volume. 
In accordance with the usual procedure, we take the window function to be a top hat. 
We normalize a{R,z) to Ug = (T(8/hMpc,0) = 0.9 ±0.1 [T], taking into account the 
suppression of the growth of A at late times due to the cosmological constant (see section 
16.4|) . as our analytic calculation of the transfer function breaks down for modes close 
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Figure 6. The processed power spectra for the three benchmark inflation models and 
also a scale Invariant primordial power spectrum for WIMP benchmark C (from top 
to bottom: hybrid inflation, scale invariant primordial spectrum, power law and 
chaotic Inflation). As before the dotted lines are the power spectra without the effects 
of collisional damping and free streaming and the vertical dashed line denotes kf,. 

to fceq- For the purpose of estimating z^i we ignore the effects of baryons (see section 
16.11 and ESD and assume a A plus matter universe with and Ha as determined by 
WMAP [T]. Specihcally, we take = 0.370. 

Figures [3 and |H1 show Zni, as dehned by equation (ED, as a function of the scale 
R, for the benchmark WIMPs and primordial power spectra respectively. In each case 
the cut-off in the processed power spectrum at A; ~ lO^Mpc”^ leads to a plateau with 
^ni = 2:^“ at i? < i?min = 0{l) pc. 

We can now give a more precise picture of the onset of the hierarchical structure 
formation process; non-linear structure formation starts at a redshift For the 

benchmark WIMP models we see that the order of magnitude variation in /cfg leads to 
a similar variation in i?min and also (because of the dependence of the amplitude of the 
peak of the power spectrum on the cut-off scale, see hgureE)) a range of values z^{^^ ~ 50 
to 65. For the scale-dependent primordial power spectra the factor of 1.6 (see section I772|l 
variation in the amplitude of the fluctuations on scales k ~ 10®Mpc“^ translates directly 
into a comparable range of values for Thus, for plausible WIMP properties and for 
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Figure 7. The redshift at which typical fluctuations of comoving scale R become 
non-linear, Zni, for the WIMP benchmark models discussed in the text (from top to 
bottom: D, C, B, A). The full lines take into account the effects of collisional damping 
and free streaming, whereas the dashed line shows the behaviour without a cut-off in 
the power spectrum. The normalisation is fixed by as = 0.9. 


a range of inflation models which produce scale-dependent power spectra consistent with 
the WMAP data for typical fluctuations takes values in the range 40 to 80. For 
the best fit WMAP matter density and a scale invariant power spectrum = 60± 10 
with the variation being due to the dependence of the free-streaming scale on the WIMP 
properties. 

The redshift of formation of the very first WIMPy halos is very different from 
the redshift z^^^^ when hierarchical structure formation starts at a typical place in the 
universe. In order to make quantitative statements about the hrst WIMPy halos we 
need to specify the statistical distribution of density fluctuations. Here we assume that 
they have a normal (gaussian) distribution, which is justihed by two physical arguments. 
Firstly, we are looking at the very hrst non-linear objects to form, thus there was no 
non-linear physics before the rare huctuations that we are going to discuss enter the 
non-linear regime. This means that, if the primordial huctuations are gaussian, the 
huctuations from which the hrst non-linear (rare) objects form should also be close to 
gaussian. On top of that argument, the central limit theorem (in the limit of large 
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Figure 8. The redshift at which typical fluctuations of comoving scale R become 
non-linear, Zni, now for the primordial power spectra discussed in the text (from top 
to bottom: hybrid inflation, scale invariant primordial spectrum, power law and rn?(jP 
chaotic inflation). The normalisation is fixed by erg = 0.9. 


numbers of distributions - one for each trial volume - the cumulative distribution 
becomes normal) justihes the use of a gaussian distribution. A typical comoving volume 
of interest is the mass collection volume of the Milky Way, which is about 1 Mpc^. As 
shown above, the cut-off due to free-streaming implies that the hrst WIMPy halos have 
a mass collecting volume of about 1 pc^. Thus we are talking about a sample of ~ 10^® 
primary halos. 

As a consequence the probability that any one of the 10^® primary regions within 
a comoving Milky Way volume has a density that exceeds Na is given by 


P(A^(x) > Na) 


Such a primary halo with 


= - l-eii{N/V2) . 

Na goes non-linear at 


( 80 ) 




max 

^nl 


m 


1 . 


(81) 


As fluctuations grow linear with the scale factor during the matter dominated epoch, this 
condition provides us with the simple result ~ = 1) ^ (60 ± 20)N. 

Here we consider the uncertainty in the primordial power spectrum and in the WIMP 
physics. The comoving cut-off scale Pmin is independent of the redshift z^i- 
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We estimate the size and mass of the hrst generation of subhalos that form at 
as well as the size and mass of the very hrst WIMPy halos, using the spherical collapse 
model (see e.g. Reference EDI). We should caution that this simplihed model has not 
yet been validated in this regime where the scale dependence of the (processed) power 
spectrum is relatively weak. The mean CDM mass within a sphere of comoving radius R 
is M{R) = 1.6 X lO“^M0(a;m/O.14)(R/pc)^. CDM overdensities that go non-linear have 
mass twice this value i.e. roughly equal to the mass of Mars. These WIMP halos are 
however much less compact than Mars. The physical size of the hrst halos at turn-around 
(when their evolution decouples from the cosmic expansion) is r = 1.05i?/[l -|- 
which is ~ (0.02/A^) pc for Rmin = 1 pc. The hrst halos then undergo violent relaxation, 
decreasing in radius by a factor of two so that their present day radius would be of order 
tens of milli-pc {N = 1), comparable to the size of the solar system, and smaller. 

A rough estimate of the relevance and chances of survival of the very hrst WIMPy 
halos can be made using the present day density contrast of these objects in the spherical 
collapse model. We hnd 
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An 
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2M(R„ 


r(i?n 


qv)) 


1 3 


CcO 


3.7(60 ± 20)^iV^ 


(82) 


Remarkably, the result is independent of the comoving size As larger halos form 

later, their density contrast is smaller, as A oc hgure IHl we plot A(iV) for 

z^^^ = 40, 60 and 80 and compare it to the density contrast in the galactic disc and in the 
halo in the solar neighborhood, Adisc = (0.3 to 1.2)10® and Ahaio = (0.2 to 1.3)10® jlHl- 
For the typical huctuations {N = 1) we hnd A = (0.2 to 1.8)10®, which is of the 
same order of magnitude as the local density contrast of the disc. This suggests that 
the typical hrst WIMPy halos may not survive the highly non-linear processing which 
occurs during structure formation. 

The situation is very diherent if we look at the rare huctuations {N > 1), which are 
characterised by the fact that they go nonlinear much earlier and are thus much denser 
than other close-by structures. As can be seen in hgure El e.g. N = 3 overdense regions 
lead to a range of A{N = 3) = (0.6 to 4.9)10^, more than an order of magnitude denser 
than the local disc and more than two orders of magnitude above the local halo density. 
Statistically it needs ~ 740 comoving pc® volumes to hnd one N = 3 huctuation. The 
Milky Way has a (comoving) mass collecting radius of ~ IMpc therefore, if these very 
hrst WIMP halos survive, there would be roughly 10®®(10®®) ‘rare N = 3(6)’ WIMPy 
subhalos within the Milky Way. Assuming the Milky Way has a volume of the order of 
(100 kpc)®, there would be, on average, one ‘rare’ N = 3(6) subhalo in each pc® [(20 
pc)®] volume. 

What are the possible implications for direct and indirect dark matter searches? 
As we discussed above, the comparison of the local halo and disc density contrasts with 
those expected for the hrst typical huctuations indicates that most of the hrst generation 
of halos will likely be destroyed during the structure formation process. In this case, 
apart from within the few rare surviving subhalos, the direct detection event rate will 
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Figure 9. The present day (mean) density contrast of Na fluctuations A(A^) for, 
from top to bottom, = 1) = 80,60,40. The dashed (dotted) lines indicate 

the range of values for the density contrast of the Milky Way disc (halo) in the solar 
neighbourhood. The normalisation is fixed by as = 0.9. 


be only slightly lower than found using the standard approach, which assumes a smooth 
dark matter distribution. The consequences for indirect detection could be much more 
dramatic. As the indirect detection rate scales with the square of the density contrast, 
the ‘rare’ subhalos discussed above could provide ‘bright’ point sources in the solar 
neighborhood, which eventually could dominate other sources (e.g. the center of the 
Milky Way). A detailed investigation of the effect on the expected direct and indirect 
dark matter search rates, as well as a detailed investigation of the survival probability 
of rare fluctuations, is beyond the scope of this paper. 

9. Discussion 

Extensive experimental efforts are being devoted to detecting WIMPs, either directly 
in the lab or indirectly via their annihilation products. Results from these searches are 
usually quoted in terms of limits on the relevant WIMP interaction cross section. The 
expected event rates in fact depend, in some cases very sensitively, on the dark matter 
distribution and hence the formation history of the Milky Way. 
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Collisional damping and free-streaming produce a cut-off in the power spectrum 
at a co-moving scale around 1 pc and set the scale of the hrst, smallest WIMP halos 
to from. In this paper we have calculated the damping processes for generic WIMPs 
and examined the effect of scale dependence of the primordial perturbation spectrum. 
We have hxed the parameters of the ACDM model, a;m,a;b, flv and A (resp. us), to the 
best-£t WMAP values. A discussion of the effect of varying CDM density has been 
provided in our previous work [Hj. 

We have found that the smallest scale fluctuations go non-linear (more precisely 
typical one-sigma fluctuations collapse to from dark matter halos) at a red-shift in the 
range 40-80, with the hrst WIMP halos forming signihcantly earlier from rare large 
huctuations. Finally we estimated the properties of both the typical and hrst small 
halos to form using the spherical collapse model. The mass of the halos is independent 
of the size of the huctuation from which they form, however the hrst rare huctuations 
to collapse form more compact, denser halos. 

Neglecting collisional damping and free streaming of WIMPs would result in 
monotonically increasing power of density huctuations on small scales. As a 
consequence, there would be a divergence of the energy density of the huctuations at 
small scales and some kind of regularization procedure would be required to make the 
hierarchical picture of structure formation well-dehned. The collisional damping and 
free streaming of WIMPs regularize the power spectrum by providing a physical cut-oh 
scale. 

Numerical studies of the formation of the hrst WIMPy halos have recently been 
carried out by Diemand et ah The subsequent evolution of these halos, and the 
resulting present day dark matter distribution remains an important outstanding issue. 
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